ADIM: Bayesian inference using the integrated nested Laplace approximation (INLA)

.title[
# ADIM: Bayesian inference using the integrated nested Laplace approximation (INLA)
]
.subtitle[
## Master’s Degree in Data Analysis, Process Improvement and Decision Support Engineering
]
.author[
### Joaquín Martínez-Minaya, 2024-12-09 <a href="http://vabar.es/" style="color:white;"> VAlencia BAyesian Research Group </a> <a href="https://smeg-bayes.org/ " style="color:white;"> Statistical Modeling Ecology Group </a> <a href="https://giem.blogs.upv.es/" style="color:white;"> Grupo de Ingeniería Estadística Multivariante </a> 
]
.date[
### <a href="jmarmin@eio.upv.es" style="color:white;"> jmarmin@eio.upv.es </a>
]

---

#Outline

## 1. Why INLA?

## 2. Elements to understand how INLA works

## 3. Putting all the pieces together: INLA

## 4. `R-INLA`

## 5. Model Selection

## 6. Examples
## 7. References

---
class: center, middle, animated, rotateInUpRight, inverse

# 1. Why INLA?

---
# INLA as an alternative to MCMC

- MCMC is an asymptotically exact method whereas INLA is an .hlb[approximation]. Their error are frequently very similar, as has been shown in many simulation studies.

- INLA is a .hlb[fast alternative] to MCMC for the general class of latent Gaussian models (LGMs). Many familiar models can be re-cast to look like LGMs:
  - .hlb[generalized linear models], .hlb[generalized additive models], smoothing spline models, 
  - state space models, semi-parametric regression, .hlb[random walk (first and second order)] models, longitudinal data models,
  - .hlb[spatial and spatiotemporal] models, log-Gaussian Cox processes and geostatistical and geoadditive models., etc.
--

- To understand INLA, we need to be familiar with:
  - Latent Gaussian models
  - Gaussian Markov Random Fields (GMRFs)
  - Laplace approximations
  
  
  
---

# 2. Elements to understand how INLA works

---

# Latent Gaussian model

## Level 1 : likelihood
The first stage is formed by the .hlb[conditionally independent likelihood] function of data coming from a certain exponential family distribution:

`$$p(\boldsymbol{y} \mid \boldsymbol{\theta}, \boldsymbol{\psi}_1) = \prod_{i=1}^{n} p(y_i \mid \eta_i(\boldsymbol{\theta}),\boldsymbol{\psi}_1)$$`
- `$\boldsymbol{y} = (y_1,\ldots,y_{n})^T$` is the response vector, `$\boldsymbol{\theta} = (\theta_1,\ldots,\theta_n)^T$` is the .hlb[latent field],
- `$\boldsymbol{\psi}_1$` is the hyperparameter vector of the exponential family distribution and
- `$\eta_i(\boldsymbol{\theta})$` is the `$i$`-th linear predictor that connects the data to the latent field.
	
Indeed each `$\eta_i$` can take a more general additive form:
  `$$\eta_i  = \beta_0 + \sum_{j=1}^{J}\beta_k x_{ij} + \sum_{k=1}^{K}f^{(k)}(z_{ik})$$`
		
---

# Latent Gaussian model

## Level 2: latent Gaussian field
- The second stage is formed by the .hlb[latent Gaussian field], where we attribute a Gaussian distribution with mean `$\boldsymbol{\mu}$` and precision matrix `$Q(\boldsymbol{\psi}_2)$` to the latent field `$\boldsymbol{\theta}$` conditioned on the hyperparameters `$\boldsymbol{\psi}_2$`, that is:
$$ \boldsymbol{\theta} \mid \boldsymbol{\psi}_2 \sim \mathcal{N}(\boldsymbol{0},Q^{-1}(\boldsymbol{\psi}_2)) $$

- If we can assume conditional independence in `$\boldsymbol{\theta}$`, then this latent field is a .hlb[Gaussian Markov Random Field (GMRF)].

## Level 3: hyperparameters
- Finally, the third stage is formed by the .hlb[prior distribution] assigned to the hyperparameters:
`$$\boldsymbol{\psi} = (\boldsymbol{\psi}_1,\boldsymbol{\psi}_2) \sim p(\boldsymbol{\psi})$$`

---
# Mixed Models Using INLA

- A study investigates agreement between devices measuring **respiratory rates** in **COPD patients** across **11 activities**, comparing a **chest-band device** to a **gold standard device** (Oxycon mobile).  
- The dataset includes 21 subjects performing activities such as sitting, walking, and climbing stairs.

The variables are:
- **y**: respiratory rate (breaths per minute).  
- **device**: measurement device (**oxicon**, **chest_band**).  
- **replicate**: replicate measurements within each activity/subject.  
- **act**: activity type (11 levels, e.g., sitting, climbing stairs).

| subj|        y|replicate |act     |device |
|----:|--------:|:---------|:-------|:------|
|    1| 38.19294|1         |Sitting |oxicon |
|    1| 40.65189|2         |Sitting |oxicon |
|    1| 36.00310|3         |Sitting |oxicon |
|    1| 32.39582|1         |Lying   |oxicon |
|    1| 24.67484|2         |Lying   |oxicon |
|    1| 26.30875|3         |Lying   |oxicon |

---
# Gaussian Markov Random Fields (GMRFs)

- A GMRF is a random vector following a .hlb[multivariate normal distribution] with Markov properties.
`$$i \neq j, \ \theta_i \mid \theta_{ij},$$`
being `$-ij$` all elements other than `$i$` and `$j$`.

- Rue et al. (2009) showed how .hbl[conditional independence properties] are encoded in the precision matrix, and how this can be exploited to improve computation involving these matrices. 
`$$i \neq j, \ \theta_i \perp \theta_j \mid \theta_{ij},$$`
`$$\theta_i \perp \theta_j \mid \boldsymbol{\theta}_{ij} \leftrightarrow \boldsymbol{Q}_{ij} = 0$$`
--

- This Markov assumption in the GMRF results in a .hlb[sparse precision matrix]. This sparseness aids extremely fast computation.

---
# The pairwise Markov property
## The two black nodes are conditionally independent given the gray nodes

---
# Example: precision matrix in AR1

|       |       |       |       |       |
|------:|------:|------:|------:|------:|
| 0.8730| 0.6957| 0.5201| 0.3460| 0.1728|
| 0.6957| 1.3931| 1.0417| 0.6929| 0.3460|
| 0.5201| 1.0417| 1.5659| 1.0417| 0.5201|
| 0.3460| 0.6929| 1.0417| 1.3931| 0.6957|
| 0.1728| 0.3460| 0.5201| 0.6957| 0.8730|
]

|        |        |        |        |        |
|-------:|-------:|-------:|-------:|-------:|
|  1.9025| -0.9500|  0.0000|  0.0000|  0.0000|
| -0.9500|  1.9025| -0.9500|  0.0000|  0.0000|
|  0.0000| -0.9500|  1.9025| -0.9500|  0.0000|
|  0.0000|  0.0000| -0.9500|  1.9025| -0.9500|
|  0.0000|  0.0000|  0.0000| -0.9500|  1.9025|
]

---
# Laplace approximations
- .hlb[The Laplace approximation] is used to estimate any distribution `$p(\theta)$` with a normal distribution.

- It uses the first three terms (quadratic function) .hlb[Taylor series expansion] around the mode `$\theta^*$` of a function to approximate its log.

- Using the approximation, `$p(\theta)$` can be approximated using a .hlb[Gaussian distribution] with mean the mode `$\theta^*$` and variance the Fisher information, `$\frac{-1}{\frac{d^2 \log(p(\theta^*))}{d \theta^2}}$`.

`$$p(\theta) \approx \mathcal{N} \left(\theta^*, \frac{-1}{\frac{d^2 \log(p(\theta^*))}{d \theta^2}} \right)$$`

- It can be easily expanded to the multivariate case.

---
# Example: approximating the beta distribution

---
class: inverse, center, middle, animated, bounceInDown

# 3. Putting all the pieces together: INLA

---
# INLA: Aim

## Marginals of the latent field and hyperparameters
`$${p(\theta_i \mid \boldsymbol{y})} =  \int {p(\theta_i \mid \boldsymbol{\psi}, \boldsymbol{y})} \cdot {p (\boldsymbol{\psi} \mid \boldsymbol{y})} d \boldsymbol{\psi}\,\,, \ i=1, \ldots, n$$`
`$$p(\psi_j \mid \boldsymbol{y}) = \int{p(\boldsymbol{\psi} \mid \boldsymbol{y})}  d \boldsymbol{\psi}_{-j}\,, \ j=1, \ldots, m$$`

- As a result, we have to numerically approximate:

1. The .hlb[joint posterior distribution of the hyperparmeters] `${p(\boldsymbol{\psi} \mid \boldsymbol{y})}$`, needed to calculate the posterior hyperparameters marginals `$p(\psi_j \mid \boldsymbol{y})$`, and the .hbl[posterior marginals of the latent field] `$p(\theta_i \mid  \boldsymbol{y})$`.
	
	2. The .hlb[marginals of the full conditional distribution] of `$\boldsymbol{\theta}$`, `${p(\theta_i \mid \boldsymbol{\psi}, \boldsymbol{y})}$`, needed to compute the posterior marginals of the latent field `$p(\theta_i \mid \boldsymbol{y})$`.

---
# Hyperparameters: joint posterior distribution
- The approximation is computed as follows

`$$\tilde{p}(\boldsymbol{\psi} \mid \boldsymbol{y}) := \frac{p(\boldsymbol{\theta},\boldsymbol{\psi}| \boldsymbol{y})}{p_G(\boldsymbol{\theta} \mid \boldsymbol{\psi}, \boldsymbol{y})}\bigg | _{\boldsymbol{\theta}=\boldsymbol{\theta}^*(\boldsymbol{\psi})}\,\,,$$`

- where:
  - `$p_G(\boldsymbol{\theta} \mid \boldsymbol{\psi}, \boldsymbol{y})$` is the Gaussian approximation to the full conditional of `$\boldsymbol{\theta}$`, `$p(\boldsymbol{\theta} \mid \boldsymbol{\psi}, \boldsymbol{y})$` given by the .hlb[Laplace method], and,
  
  - `$\boldsymbol{\theta}^*(\boldsymbol{\psi})$` is the mode of the full conditional of `$\boldsymbol{\theta}$` for a given `$\boldsymbol{\psi}$`.

- Note: this approximation is exact if `$p(\boldsymbol{\theta}\mid \boldsymbol{y}, \boldsymbol{\psi})$` is Gaussian.

---
# Full posterior marginals for the latent field

### .hlb[Gaussian approximation] 
- Conditional posterior distributions `${p(\theta_i \mid \boldsymbol{\psi}, \boldsymbol{y})}$` are approximated directly as the marginals from `$p_G(\boldsymbol{\theta} \mid \boldsymbol{\psi}, \boldsymbol{y})$`.
- It is the .hlbred[fastest to compute] but with possible .hlbred[errors] in the location of the posterior mean.

--
	
### .hlb[Laplace approximation] 
- The vector `$\boldsymbol{\theta}$` is rewriten as `$\boldsymbol{\theta}=(\theta_i, \boldsymbol{\theta}_{-i})$`, and the Laplace approximation is used for each element of the latent field 
	`$$\tilde{p}(\theta_i \mid \boldsymbol{\psi}, \boldsymbol{y}) := \frac{p(\boldsymbol{\theta},\boldsymbol{\psi}| \boldsymbol{y})}{p_{LG}(\boldsymbol{\theta}_{-i} \mid \theta_i ,\boldsymbol{\psi}, \boldsymbol{y})}\bigg |_{\boldsymbol{\theta}_{-i}=\boldsymbol{\theta}_{-i}^*(\theta_i, \boldsymbol{\psi})}\,\,,$$`
	where `$p_{LG}(\boldsymbol{\theta}_{-i} \mid \theta_i ,\boldsymbol{\psi}, \boldsymbol{y})$` is the Laplace Gaussian approximation to `$p(\boldsymbol{\theta}_{-i} \mid \theta_i ,\boldsymbol{\psi}, \boldsymbol{y})$` and `$\boldsymbol{\theta}_{-i}$` is its mode.
	
- The .hlbred[most accurate] but .hlbred[time consuming].

---

# Full posterior marginals for the latent field
### .hlb[Simplified Laplace approximation]

- Based on a Taylor's series expansion of third order.
- .hlbred[Fast to compute] and usually .hlbred[accurate enough].

---

# Final step: integration
- The INLA algorithm uses Newton-like methods to explore the joint posterior distribution for the hyperparameters `$\tilde{p}(\boldsymbol{\psi}|\boldsymbol{y})$` to find .hlb[suitable points] for the numerical integration.

- Posterior marginals for the .hlb[latent variables] `$\tilde{p}(\theta_i|\boldsymbol{y})$` are then computed via numerical integration as: 
  `$$\tilde{p}(\theta_i \mid \boldsymbol{y}) = \int \tilde{p}(\theta_i \mid \boldsymbol{\psi},\boldsymbol{y}) \tilde{p}(\boldsymbol{\psi} \mid \boldsymbol{y}) \text{d}\boldsymbol{\psi} \approx \sum_{k=1}^{K} \tilde{p}(\theta_i \mid \boldsymbol{\psi}^{(k)},\boldsymbol{y}) \tilde{p}(\boldsymbol{\psi}^{(k)} \mid \boldsymbol{y})\Delta_k$$`

- Posterior marginals for the .hlb[hyperparameters] `$\psi_j$` are approximated using the integrations points previously constructed.

---

# The new avenue in INLA

- Posterior means of `$\boldsymbol{\theta}$` and `$\boldsymbol{\eta}$` might be inaccurate based on the Gaussian assumption of the conditional posterior.

- .hlb[Variational Bayes correction] to Gaussian means by [Van Niekerk and Rue (2021)](https://arxiv.org/pdf/2111.12945.pdf) can be used to eﬃciently calculate an improved mean for the marginal posteriors of the linear predictors, by improving the posterior means of the latent ﬁeld.

- It is improved based on the following variational function:

`$$E_{q(\boldsymbol{\theta} \mid \boldsymbol{y})}[-\text{log}(p(\boldsymbol{\theta}\mid \boldsymbol{y}))] + KLD[q(\boldsymbol{\theta} \mid \boldsymbol{y}) \mid \mid p(\boldsymbol{\theta})]] \,,$$`

1. `$q(.)$` a member of the variational class,

2. `$p(\boldsymbol{\theta} \mid \boldsymbol{y})$` the posterior distribution of the latent field, and

3. `$p(\boldsymbol{\theta})$` the posterior distribution of the latent field.

---
class: inverse, center, middle, animated, bounceInDown

# 4. R-INLA

---
# Data
.left-column2[

|        y|         x| id|
|--------:|---------:|--:|
| 2.109177| 0.3077661|  1|
| 1.954976| 0.2576725|  2|
| 2.294048| 0.5523224|  3|
| 2.217938| 0.0563832|  4|
| 2.007082| 0.4685493|  5|
| 2.339932| 0.4837707|  6|
]

.right-column2[
![](data:image/png;base64,#S2_2_inla_files/figure-html/unnamed-chunk-7-1.png)
]

---
# Fitting the model using `R-INLA`

## Defining the formula

``` r
formula <- y ~ 1 + x # 1 is refered to the intercept term
formula <- y ~ 1 + f(x, model = "linear")
```

## Calling `R-INLA`

``` r
model1 <- inla(formula, 
 family = 'gaussian', 
 data = data,
 control.inla = list(strategy = 'simplified.laplace'))
```
---
# Posterior distributions

## .hlb[Posterior distribution of the parameters]

|            |   mean|     sd| 0.025quant| 0.5quant| 0.975quant|   mode| kld|
|:-----------|------:|------:|----------:|--------:|----------:|------:|---:|
|(Intercept) | 1.9946| 0.0225|     1.9503|   1.9946|     2.0389| 1.9946|   0|
|x           | 0.4935| 0.0388|     0.4174|   0.4935|     0.5697| 0.4935|   0|

## .hlb[Posterior distributions of the hyperparameters]

|                                        |    mean|      sd| 0.025quant| 0.5quant| 0.975quant|    mode|
|:---------------------------------------|-------:|-------:|----------:|--------:|----------:|-------:|
|Precision for the Gaussian observations | 99.4863| 14.0682|    73.8685|  98.8272|   128.9193| 97.5018|

---
# Families

``` r
  inla.list.models(section = "likelihood")
```

```
## Section [likelihood]
##      0binomial                     New 0-inflated Binomial                 
##      0binomialS                    New 0-inflated Binomial Swap            
##      0poisson                      New 0-inflated Poisson                  
##      0poissonS                     New 0-inflated Poisson Swap             
##      agaussian                     The aggregated Gaussian likelihoood     
##      bell                          The Bell likelihood                     
##      beta                          The Beta likelihood                     
##      betabinomial                  The Beta-Binomial likelihood            
##      betabinomialna                The Beta-Binomial Normal approximation likelihood
##      bgev                          The blended Generalized Extreme Value likelihood
##      binomial                      The Binomial likelihood                 
##      cbinomial                     The clustered Binomial likelihood       
##      cennbinomial2                 The CenNegBinomial2 likelihood (similar to cenpoisson2)
##      cenpoisson                    Then censored Poisson likelihood        
##      cenpoisson2                   Then censored Poisson likelihood (version 2)
##      circularnormal                The circular Gaussian likelihoood       
##      coxph                         Cox-proportional hazard likelihood      
##      dgp                           Discrete generalized Pareto likelihood  
##      exponential                   The Exponential likelihood              
##      exponentialsurv               The Exponential likelihood (survival)   
##      fl                            The fl likelihood                       
##      fmri                          fmri distribution (special nc-chi)      
##      fmrisurv                      fmri distribution (special nc-chi)      
##      gamma                         The Gamma likelihood                    
##      gammacount                    A Gamma generalisation of the Poisson likelihood
##      gammajw                       A special case of the Gamma likelihood  
##      gammajwsurv                   A special case of the Gamma likelihood (survival)
##      gammasurv                     The Gamma likelihood (survival)         
##      gaussian                      The Gaussian likelihoood                
##      gaussianjw                    The GaussianJW likelihoood              
##      gev                           The Generalized Extreme Value likelihood
##      ggaussian                     Generalized Gaussian                    
##      ggaussianS                    Generalized GaussianS                   
##      gompertz                      gompertz distribution                   
##      gompertzsurv                  gompertz distribution                   
##      gp                            Generalized Pareto likelihood           
##      gpoisson                      The generalized Poisson likelihood      
##      iidgamma                      (experimental)                          
##      iidlogitbeta                  (experimental)                          
##      loggammafrailty               (experimental)                          
##      logistic                      The Logistic likelihoood                
##      loglogistic                   The loglogistic likelihood              
##      loglogisticsurv               The loglogistic likelihood (survival)   
##      lognormal                     The log-Normal likelihood               
##      lognormalsurv                 The log-Normal likelihood (survival)    
##      logperiodogram                Likelihood for the log-periodogram      
##      mgamma                        The modal Gamma likelihood              
##      mgammasurv                    The modal Gamma likelihood (survival)   
##      nbinomial                     The negBinomial likelihood              
##      nbinomial2                    The negBinomial2 likelihood             
##      nmix                          Binomial-Poisson mixture                
##      nmixnb                        NegBinomial-Poisson mixture             
##      npoisson                      The Normal approximation to the Poisson likelihood
##      nzpoisson                     The nzPoisson likelihood                
##      poisson                       The Poisson likelihood                  
##      poisson.special1              The Poisson.special1 likelihood         
##      pom                           Likelihood for the proportional odds model
##      qkumar                        A quantile version of the Kumar likelihood
##      qloglogistic                  A quantile loglogistic likelihood       
##      qloglogisticsurv              A quantile loglogistic likelihood (survival)
##      rcpoisson                     Randomly censored Poisson               
##      simplex                       The simplex likelihood                  
##      sn                            The Skew-Normal likelihoood             
##      stdgaussian                   The stdGaussian likelihoood             
##      stochvol                      The Gaussian stochvol likelihood        
##      stochvolln                    The Log-Normal stochvol likelihood      
##      stochvolnig                   The Normal inverse Gaussian stochvol likelihood
##      stochvolsn                    The SkewNormal stochvol likelihood      
##      stochvolt                     The Student-t stochvol likelihood       
##      t                             Student-t likelihood                    
##      tpoisson                      Thinned Poisson                         
##      tstrata                       A stratified version of the Student-t likelihood
##      tweedie                       Tweedie distribution                    
##      weibull                       The Weibull likelihood                  
##      weibullsurv                   The Weibull likelihood (survival)       
##      wrappedcauchy                 The wrapped Cauchy likelihoood          
##      xbinomial                     The Binomial likelihood (expert version)
##      xpoisson                      The Poisson likelihood (expert version) 
##      zeroinflatedbetabinomial0     Zero-inflated Beta-Binomial, type 0     
##      zeroinflatedbetabinomial1     Zero-inflated Beta-Binomial, type 1     
##      zeroinflatedbetabinomial2     Zero inflated Beta-Binomial, type 2     
##      zeroinflatedbinomial0         Zero-inflated Binomial, type 0          
##      zeroinflatedbinomial1         Zero-inflated Binomial, type 1          
##      zeroinflatedbinomial2         Zero-inflated Binomial, type 2          
##      zeroinflatedcenpoisson0       Zero-inflated censored Poisson, type 0  
##      zeroinflatedcenpoisson1       Zero-inflated censored Poisson, type 1  
##      zeroinflatednbinomial0        Zero inflated negBinomial, type 0       
##      zeroinflatednbinomial1        Zero inflated negBinomial, type 1       
##      zeroinflatednbinomial1strata2 Zero inflated negBinomial, type 1, strata 2
##      zeroinflatednbinomial1strata3 Zero inflated negBinomial, type 1, strata 3
##      zeroinflatednbinomial2        Zero inflated negBinomial, type 2       
##      zeroinflatedpoisson0          Zero-inflated Poisson, type 0           
##      zeroinflatedpoisson1          Zero-inflated Poisson, type 1           
##      zeroinflatedpoisson2          Zero-inflated Poisson, type 2           
##      zeroninflatedbinomial2        Zero and N inflated binomial, type 2    
##      zeroninflatedbinomial3        Zero and N inflated binomial, type 3
```

---
# Latent effects

``` r
  inla.list.models(section = "latent")
```

```
## Section [latent]
##      2diid                         (This model is obsolute)                
##      ar                            Auto-regressive model of order p (AR(p))
##      ar1                           Auto-regressive model of order 1 (AR(1))
##      ar1c                          Auto-regressive model of order 1 w/covariates
##      besag                         The Besag area model (CAR-model)        
##      besag2                        The shared Besag model                  
##      besagproper                   A proper version of the Besag model     
##      besagproper2                  An alternative proper version of the Besag model
##      bym                           The BYM-model (Besag-York-Mollier model)
##      bym2                          The BYM-model with the PC priors        
##      cgeneric                      Generic latent model specified using C  
##      clinear                       Constrained linear effect               
##      copy                          Create a copy of a model component      
##      crw2                          Exact solution to the random walk of order 2
##      dmatern                       Dense Matern field                      
##      fgn                           Fractional Gaussian noise model         
##      fgn2                          Fractional Gaussian noise model (alt 2) 
##      generic                       A generic model                         
##      generic0                      A generic model (type 0)                
##      generic1                      A generic model (type 1)                
##      generic2                      A generic model (type 2)                
##      generic3                      A generic model (type 3)                
##      iid                           Gaussian random effects in dim=1        
##      iid1d                         Gaussian random effect in dim=1 with Wishart prior
##      iid2d                         Gaussian random effect in dim=2 with Wishart prior
##      iid3d                         Gaussian random effect in dim=3 with Wishart prior
##      iid4d                         Gaussian random effect in dim=4 with Wishart prior
##      iid5d                         Gaussian random effect in dim=5 with Wishart prior
##      iidkd                         Gaussian random effect in dim=k with Wishart prior
##      intslope                      Intecept-slope model with Wishart-prior 
##      linear                        Alternative interface to an fixed effect
##      log1exp                       A nonlinear model of a covariate        
##      logdist                       A nonlinear model of a covariate        
##      matern2d                      Matern covariance function on a regular grid
##      meb                           Berkson measurement error model         
##      mec                           Classical measurement error model       
##      ou                            The Ornstein-Uhlenbeck process          
##      revsigm                       Reverse sigmoidal effect of a covariate 
##      rgeneric                      Generic latent model specified using R  
##      rw1                           Random walk of order 1                  
##      rw2                           Random walk of order 2                  
##      rw2d                          Thin-plate spline model                 
##      rw2diid                       Thin-plate spline with iid noise        
##      scopy                         Create a scopy of a model component     
##      seasonal                      Seasonal model for time series          
##      sigm                          Sigmoidal effect of a covariate         
##      slm                           Spatial lag model                       
##      spde                          A SPDE model                            
##      spde2                         A SPDE2 model                           
##      spde3                         A SPDE3 model                           
##      z                             The z-model in a classical mixed model formulation
```

---

# Hyperpriors

``` r
  inla.list.models(section = "prior")
```

```
## Section [prior]
## betacorrelation Beta prior for the correlation 
## dirichlet Dirichlet prior 
## expression: A generic prior defined using expressions
## flat A constant prior 
## gamma Gamma prior 
## gaussian Gaussian prior 
## invalid Void prior 
## jeffreystdf Jeffreys prior for the doc 
## laplace Laplace prior 
## linksnintercept Skew-normal-link intercept-prior 
## logflat A constant prior for log(theta) 
## loggamma Log-Gamma prior 
## logiflat A constant prior for log(1/theta) 
## logitbeta Logit prior for a probability 
## logtgaussian Truncated Gaussian prior 
## logtnormal Truncated Normal prior 
## minuslogsqrtruncnormal (obsolete) 
## mvnorm A multivariate Normal prior 
## none No prior 
## normal Normal prior 
## pc Generic PC prior 
## pc.alphaw PC prior for alpha in Weibull 
## pc.ar PC prior for the AR(p) model 
## pc.cor0 PC prior correlation, basemodel cor=0 
## pc.cor1 PC prior correlation, basemodel cor=1 
## pc.dof PC prior for log(dof-2) 
## pc.fgnh PC prior for the Hurst parameter in FGN 
## pc.gamma PC prior for a Gamma parameter 
## pc.gammacount PC prior for the GammaCount likelihood 
## pc.gevtail PC prior for the tail in the GEV likelihood
## pc.matern PC prior for the Matern SPDE 
## pc.mgamma PC prior for a Gamma parameter 
## pc.prec PC prior for log(precision) 
## pc.range PC prior for the range in the Matern SPDE
## pc.sn PC prior for the skew-normal 
## pc.spde.GA (experimental) 
## pom #classes-dependent prior for the POM model
## ref.ar Reference prior for the AR(p) model, p<=3
## rprior: A R-function defining the prior 
## table: A generic tabulated prior 
## wishart1d Wishart prior dim=1 
## wishart2d Wishart prior dim=2 
## wishart3d Wishart prior dim=3 
## wishart4d Wishart prior dim=4 
## wishart5d Wishart prior dim=5 
## wishartkd Wishart prior
```

---

# 5. Model Selection

---
# Model selection scores in `R-INLA`
- When use different covariates and random effects, we need some measures to select the best model:
  - .hlb[DIC]: deviance information criteria
		`$$DIC  = 2*E(D(\boldsymbol{\theta})) - 						D(E(\boldsymbol{\theta}))$$`
		
	- .hlb[WAIC]: within-sample predictive score
		`$$WAIC  = \sum_{i} var_{post}(log( p(y_i|\boldsymbol{\theta})))$$`
		
	- .hlb[LCPO]: leave-one-out cross-validation score
	  `$$CPO_i = p(y_i \mid y_{-i})$$`
		
	  `$$LCPO  = - \overline{\log(CPO_i) }$$`

---
class: inverse, center, middle, animated, rotateInDownRight

# 6. Example

---
# Bayesian splines
- GLMM with independent random effect does not cover situations in which relationship between the response variable and the covariate is not linear.

- In INLA, we can do this by means of the .hlb[random walk] of order 1 and 2.
  - .hlb[First order Random Walk (RW1)] 
`$$\Delta x_j = x_j - x_{j+1} \sim  \mathcal{N}\left(0, \sigma^2=\frac{1}{\tau}\right)$$`

- .hlb[Second order Random Walk (RW2)] `$$\Delta^2x_i = x_{i} - 2x_{i+1} + x_{i+2} \sim  \mathcal{N}\left(0, \sigma^2=\frac{1}{\tau}\right)$$`
  - The prior for the hyperparameter `$\tau$` is reparametrized in terms of their logarithm:
		`$$\log(\tau) \sim \text{logGamma}(1, 5 \cdot 10^{-5})\,\,.$$`

---

# Smoothing time series of binomial data
.left-column3[
- The number of .hlb[occurrences of rainfall] over 1 mm in the Tokyo area for each calendar year during two years (1983-84) are registered.

- It is of interest to estimate the underlying probability `${\pi_t}$` of rainfall for calendar day `${t}$` which is, a priori, assumed to change gradually over time.

- For each day `${t = 1, ..., 366}$` of the year we have the number of days that rained `${y_t}$` and the number of days that were observed `${n_t}$`.
]

**Dataset**

|  y|  n| time|
|--:|--:|----:|
|  0|  2|    1|
|  0|  2|    2|
|  1|  2|    3|
|  1|  2|    4|
|  0|  2|    5|
|  1|  2|    6|

]
---
# Smoothing time series of binomial data. The model
- A conditionally independent .hlb[binomial likelihood] function:
    $$y_t \mid \pi_t \sim \text{Binomial}(n, \pi_t),\ t = 1, ..., 366 $$
with (usual) logit link function:
    `$$\pi_t = \frac{\exp(\eta_t)}{1 + \exp(\eta_t)}$$`

- We assume that (instead of a linear predictor), `$\eta_t = f_t$`, where `$f_t$` follows a circular .hlb[random walk] of second order (RW2) model with precision `$\tau$`:
    `$$\Delta^2f_i = f_i - 2 f_{i+1} + f_{i+2} \sim \mathcal{N}(0,\tau^{-1}).$$`
The fact that we use a circular model here means that in this case `$f_1$` is a neighbor of  `$f_{366}$`. So, in this case `$\boldsymbol{\theta} = (f_1, ..., f_{366})$` and again `$\boldsymbol{\theta}|\boldsymbol{\psi}$` is .hlb[Gaussian distributed].

- To assing the prior of `$\boldsymbol{\psi} = (\tau)$`:
    `$$\log(\tau) \sim \text{logGamma}(1, 5 \cdot 10^{-5})\,\,.$$`
		
---
# Posterior distribution of the probability

---
class: inverse, center, middle, animated, slideInLeft

# 7. References

---
# This material has been constructed based on:

- Gómez-Rubio, V. (2020). Bayesian inference with INLA. Chapman and Hall/CRC.

- Parker, R. A., Scott, C., Inácio, V., & Stevens, N. T. (2020). Using multiple agreement methods for continuous repeated measures data: a tutorial for practitioners. BMC Medical Research Methodology, 20, 1-14.

- Rue, H., & Held, L. (2005). Gaussian Markov random fields: theory and applications. Chapman and Hall/CRC.

- Rue, H., Martino, S.,  & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the royal statistical society: Series b (statistical methodology), 71(2), 319-392.

- Wang, X., Ryan, Y. Y., & Faraway, J. J. (2018). Bayesian Regression Modeling with INLA. Chapman and Hall/CRC.

- <a href="https://haakonbakka.bitbucket.io/organisedtopics.html/" style="color:#FF0000;"> Tutorials by Haakon Bakka </a>

- <a href="https://www.precision-analytics.ca/blog/a-gentle-inla-tutorial/" style="color:#FF0000;"> A gentle INLA tutorial by Kathryn Morrison </a>

- <a href="https://becarioprecario.bitbucket.io/inla-gitbook/index.html" style="color:#FF0000;"> INLA book by Virgilio Gómez-Rúbio </a>

---
class: inverse, left, middle, animated, bounceInDown

# ADIM: Bayesian inference using the integrated nested Laplace approximation (INLA)
## Master's Degree in Data Analysis, Process Improvement and Decision Support Engineering

Joaquín Martínez-Minaya, 2024-12-09 
 
<a href="http://vabar.es/" style="color:white;" > VAlencia BAyesian Research Group </a> <a href="https://smeg-bayes.org/ " style="color:white;"> Statistical Modeling Ecology Group </a> <a href="https://giem.blogs.upv.es/" style="color:white;"> Grupo de Ingeniería Estadística Multivariante </a>

<a href="jmarmin@eio.upv.es" style="color:white;"> jmarmin@eio.upv.es </a>